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Using two different numerical metliods, we study the behavior of two-component Fermi gases 
interacting through short-range s-wave interactions in a harmonic trap. A correlated Gaussian basis- 
set expansion technique is used to determine the energies and structural properties, i.e., the radial 
one-body densities and pair distribution functions, for small systems with either even or odd A'^, as 
functions of the s-wave scattering length and the mass ratio k of the two species. Particular emphasis 
is put on a discussion of the angular momentum of the system in the BEC-BCS crossover regime. At 
unitarity, the excitation spectrum of the four-particle system with total angular momentum L = 
is calculated as a function of the mass ratio k. The results are analyzed from a hyperspherical 
perspective, which offers new insights into the problem. Additionally, fixed-node diffusion Monte 
Carlo calculations are performed for equal-mass Fermi gases with up to A = 30 atoms. We focus 
on the odd-even oscillations of the ground state energy of the equal-mass unitary system having up 
to A = 30 particles, which are related to the excitation gap of the system. Furthermore, we present 
a detailed analysis of the structural properties of these systems. 

PACS numbers: 



I. INTRODUCTION 

Pure Fermi systems with essentially any interaction 
strength can be realized experimentally with ultracold 
atomic gases. In most experiments to date, large sam- 
ples of atomic Li or K are trapped optically in two differ- 
ent hyperfine states, in the following simply referred to 
as "spin-up" and "spin-down" states. By tuning an ex- 
ternal magnetic field in the vicinity of a Fano-Feshbach 
resonance [l|, 0) S 3 > ^^"^ interspecies s-wave scattering 
length can be varied from non-interacting to infinitely 
strongly- interacting (either attractive or repulsive). This 
tunability is unique to atomic systems, and it has enabled 
for the first time quantitative experimental studies of the 
crossover from the molecular BEC-regime to the atomic 
BCS-regime H, i, 0, i, S ■ Since the systems studied 
experimentally are in general large, many observations 
have been explained quite successfully by applying the- 
oretical treatments based on the local density approxi- 
mation (LDA); see, e.g., Ref. [ll| and references therein. 
The LDA uses the equation of state of the homogeneous 
system as input, and, in general, accurately describes the 
properties of the system near the trap center, where the 
density changes slowly. However, it fails to accurately 
describe the properties of the system near the edge of 
the cloud, where the density varies more rapidly. 

In a different set of experiments, atomic Fermi gases 
are loaded into an optical lattice with variable barrier 
height [l^ . [isl . [T3| • In the regime where the tunneling of 
atoms between neighboring lattice sites can be neglected, 
each lattice site provides an approximately harmonic con- 
fining potential for the atoms at that site. Through the 
application of a so-called "purification scheme" [l5|, ex- 



perimentalists are now able to realize systems with a de- 
terministic number of atoms per site. So far, optical lat- 
tices have been prepared with one or zero atoms per site, 
with two or zero atoms per site, and with three or zero 
atoms per site. Optical-lattice experiments thus allow 
for the simultaneous preparation of multiple copies of 
identical few-particle systems. We anticipate that these 
experiments will be extended to larger atom samples in 
the future, thereby opening the possibility to study sys- 
tematically how the properties of the system change as 
functions of the number of atoms. Transitions from few- 
to many-body systems have, e. g., been studied experi- 
mentally in metal and rare gas clusters [H, [13 , and it is 
exciting that the experimental study of this transition in 
dilute gaseous systems is within reach. A mature body 
of theoretical work has also investigated the manner in 
which bulk electronic, magnetic and superfluid proper- 
ties can be understood by studying small or modest-size 
clusters [ll,[i3. 

This paper presents theoretical results for trapped two- 
component Fermi gases with up to A'^ = 30 fermions, 
which shed light on the few- to many-body transition 
from a microscopic or few-body point of view. To solve 
the many-body Schrodinger equation we use two differ- 
ent numerical methods, a correlated Gaussian (CG) basis 
set expansion approach and a fixed-node diffusion Monte 
Carlo (FN-DMC) approach. The CG approach allows 
for the determination of the entire energy spectrum and 
eigenstates with controlled accuracy (i.e., no approxima- 
tions are employed and the convergence can be systemat- 
ically improved) . If we demand an accuracy of the order 
of 2% or better, our current CG implementation limits 
us to treating systems with up to A^ = 6 atoms (and 
to the lowest 10 or twenty eigen states). To the best of 
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our knowledge, no other such calculations exist for dilute 
fermionic few-body systems (TV = 4 — 6) with short-range 
interactions. The FN-DMC method, in contrast, can be 
applied to larger systems but its accuracy crucially de- 
pends on the quality of the many-body nodal surface, 
which is in general unknown. Moreover, the FN-DMC 
approach as implemented here treats only ground state 
properties for the chosen symmetry. Careful comparisons 
of the ground state energy and structural properties cal- 
culated by the FN-DMC and CG approach for differ- 
ent interaction strengths validate the construction of the 
nodal surfaces employed for < 6. We expect, and pro- 
vide some evidence, that our nodal surfaces constructed 
to describe the energetically lowest-lying gas-like state of 
larger N are also quite accurate. 

Specifically, we calculate the energy of the energetically 
lowest-lying gas-like state of trapped two-species Fermi 
gases as a function of the number of particles N, the 
s-wave scattering length and the mass ratio k. Our 
ground state energies for even and odd N can be read- 
ily combined to determine the excitation gap, which is 
related to pairing physics. For small systems, we addi- 
tionally determine and discuss the excitation spectrum. 
Furthermore, we present pair correlation functions, which 
provide further insights into the pair formation process, 
and radial density profiles for the ground state. Finally, 
we elaborate on the interpretation of the behaviors within 
a framework that uses hyperspherical coordinates. This 
connection has been summarized in an earlier paper [20| . 
Here, we present additional results and discuss in more 
detail how the even-odd oscillations emerge in the hyper- 
spherical framework. Our analysis provides an alterna- 
tive means, complementary to conventional many-body 
theory, for understanding the excitation gap at unitarity. 

The remainder of this paper is organized as follows. 
Section |TT] introduces the Hamiltonian of the system un- 
der study, reviews the definitions of the normalized en- 
ergy crossover curve and the excitation gap, and sum- 
marizes some peculiar properties of the unitary gas us- 
ing hyperspherical coordinates. Section IlIII summarizes 
the CG and FN-DMC approaches, and provides some 
implementation details specific to the problem at hand. 
Section HVl presents our results for the ground state ener- 
gies, the excitation spectrum and structural properties. 
Finally, Sec. IVl concludes. 



II. THEORETICAL BACKGROUND 
A. Hamiltonian 

The main objective of this article is to obtain and 
interpret solutions to the many-body time-independent 
Schrodinger equation for a trapped two-component Fermi 
gas with short-range interactions. The model Hamilto- 
nian for iVi fermions of mass mi and A'2 fermions of mass 
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Here, fi and fii denote the position vector of the ith 
mass mi fermion and the i'th mass m2 fermion, respec- 
tively. Both atom species experience a trapping poten- 
tial characterized by the same angular frequency uj. For 
equal masses, this is indeed the case in ongoing exper- 
iments. For unequal masses, however, the two atomic 
species typically experience different trapping frequen- 
cies. Our restriction to equal trapping frequencies re- 
duces the parameter space which otherwise would be im- 
practical to explore numerically. Furthermore, our CG 
calculations simplify for equal trapping frequencies be- 
cause the center-of-mass and relative motions decouple in 
this case. The studies presented here for unequal masses 
but equal frequencies complement our earlier study 2ij, 
which treats two-component Fermi gases with unequal 
masses that experience trapping frequencies cui and uj2 
adjusted so that miWi = m2W2. In Eq. ([T]), Vb is a short- 
range two-body potential between each pair of mass mi 
and mass m2 atoms. We characterize the strength of 
Vq by the s-wave scattering length a^, which can be var- 
ied experimentally through the application of an external 
magnetic field in the vicinity of a Fano-Fcshbach reso- 
nance. Here, we model this situation by changing the 
depth of Vb; our results should be applicable to systems 
with a broad s-wave Fano-Feshbach resonance and van- 
ishingly small p-wave interactions. 

The present study considers two-component Fermi 
gases with either even or odd N, where = A^i + 
iV2. Because odd-even oscillations serve as one major sub- 
ject of this study, we set Ni = N2 for even N, and 
A^i — N2 ± 1 for odd N. In addition to the scattering 
length as, we vary the mass ratio k, 

K = mi/m2- (2) 

Throughout, we take mi > m2 so that k > 1. In most 
cases, we measure lengths in units of the oscillator length 
fl/io, o.ho ~ ■\/ft/(2/iCi;), which is defined in terms of the 
reduced mass /i, = mim2/ (mi -I- m2). 

It has been shown previously [l^, [13, [H, [l^l that small 
equal-mass two-component Fermi gases, which interact 
through short-range two-body potentials with infinitely 
large that support no s-wave bound state, support 
no tightly-bound many-body states with negative en- 
ergy. For unequal mass systems the situation is differ- 
ent [l^, i25„ J6., J7i] . Trimers consisting of two heavy par- 
ticles and one light particle that interact through short- 
range potentials support tightly-bound states with neg- 
ative energy if the mass ratio and the scattering length 
are sufficiently large. Reference [2T] discussed the role of 
non-universal trimer states for unequal-mass systems in 



3 



some detail, and we return to this discussion in Sec. lIV Al 
Throughout this work, we restrict our analysis to gas-hke 
states, consisting of atomic fermions, molecular bosons or 
both. 

To solve the Schrodinger equation for eigenstates of 
H, we use two different numerical methods: a correlated 
Gaussian (CG) basis set expansion technique and a fixed- 
node diffusion Monte Garlo (FN-DMC) technique. For 
numerical convenience, we utilize different short-range 
potentials Vq in our CG and FN-DMC calculations. We 
adopt a purely attractive Gaussian interaction potential 
defined as 



(3) 



Vb(r) = -dexp ( 



in the CG calculations, and a square well interaction po- 
tential defined as 



Voir) 



—d for r < Rq 
for r > _Ro 
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in the FN-DMC calculations. For a fixed range Rq, the 
potential depth d is adjusted so that the s-wave scatter- 
ing length Os takes the desired value. The range Rq is 
selected so that Rq <C aho- The premise is that the prop- 
erties of two-component Fermi gases with short-range in- 
teractions (or at least the universal state properties) are 
determined by the s-wave scattering length alone, and 
independent of the details of the underlying two-body 
potential if the range Rq is chosen sufficiently small. Ide- 
ally, we would consider the limit Rq — 0. This is, how- 
ever, impossible within the numerical frameworks em- 
ployed. Thus, we perform calculations for different finite 
Rq, which allows us to approximately extrapolate to the 
Rq = limit and to estimate the dependence of our re- 
sults on Rq, i.e., to estimate the scale of the finite-range 
effects. 



B. Energy crossover curve and excitation gap 

The energetically lowest-lying gas-like states of two- 
component Fermi gases with short-range interactions de- 
termine the normalized energy crossover curve and 
the excitation gap A(A^). To simplify the notation, the 
energetically lowest-lying gas-like state is referred to as 
the ground state in this section. The BCS and BEG lim- 
its of the crossover can be treated perturbatively. For 
small \as\ and as < 0, the system behaves like a weakly- 
interacting atomic Fermi gas whose leading order prop- 
erties beyond the non-interacting degenerate Fermi gas 
are determined by a^. For attractive two-body poten- 
tials that generate small as and Os > 0, in contrast, the 
system behaves like a weakly-interacting molecular Bose 
gas whose properties are to leading order determined by 
add, where add denotes the dimer-dimer scattering length. 
(One can also have small, positive with purely repul- 
sive two-body potentials that have no bound molecular 



states, but these systems behave quite differently and 
will not be considered in this paper.) In the strongly- 
interacting regime (large |as|), perturbation theory can- 
not be applied and it is not clear a priori whether the 
system behaves more like an atomic gas or a molecular 
gas, or like neither of the two. 

The definition of the normalized energy crossover curve 
introduced in Refs. [U, [2^ for even N can be 
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extended to odd N, 
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Here, E{Ni, N2) denotes the ground state energy of the 
trapped two-component gas consisting of A^i fermions 
with mass mi and N2 fermions with mass 7712- In Eq. ([5]), 
Nd is defined by 



Nd = min{iVi,iV2}, 



(6) 



and corresponds to the number of dimers formed on the 
BEG side, i.e., in the regime where as is small and posi- 
tive. Nf is defined by 



(7) 



it represents the number of unpaired atoms on the BEG 
side, and takes the value for even N and 1 for odd TV. 

In Eq. (O, En J denotes the ground state energy of the 
non-interacting two-component Fermi gas consisting of N 
atoms, where — as before — N — Ni + N2. The E^i can 
be evaluated as the sum of the noninteracting energies 
of polarized Fermi gases E^j with iVi and iV2 particles, 
Eni{N) = Elj{Ni) + El^^{N2). Following Ref. 2^9], the 
E^j{Ni) can be written in terms of the shell number n^, 
the energy of the closed shell subsystem E'^j{ns), and 
the corresponding magic number N'^^ , 



El,{K,)^E'§,{ns) + 



+ ns]{N - N''')njjj. (8) 



The shell number Us represents the number of closed 
shells and is given by 



n„ = Int 



1 



- 1 



where 



3 27A^, - J3(2437V2 _ 1; 



(9) 



(10) 



and Int[a;] is the integer part of x. Finally, the energy of 
the closed shell subsystem Ef^jius) and the correspond- 
ing magic number N'^'^ are 



j^cs^ns{ns + l){ns + 2) 
6 

EniM _ [us - l)ns{ns + l)(n. + 2) 37V^ 



(11) 
(12) 
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On the positive as side where a high-lying two-body 
bound state exists, a significant fraction of the ground 
state energy of the N fermion system is determined by 
the binding energy of the trapped dimer, which depends 

on i?o- To reduce the dependence of A^^-* ^ on the range 
i?o, the energy £'(1, 1) of Nd trapped dimer pairs is sub- 
tracted in Eq. ([5]). Thus, Aj^^^ depends to a good 
approximation only on a^, k, Ni and and not on 
the details of the underlying two-body potential (see also 
Sec. lIV A| . By construction, A^^'' changes from one on 
the weakly-interacting BCS side (small \as\ and < 0) 
to zero on the weakly-interacting molecular EEC side 
(small, positive as). 

The weakly-interacting regimes, where \ag \ <^ afi^, can 
be treated perturbatively assuming zero-range interac- 
tions, i.e., a Fermi pseudopotential [30i|. For small \as\ 
and Os < 0, the energy within first order perturbation 
theory becomes 

E^ENi + huC^^j^^^, (13) 
a-ho 

where C^^ jy, ^ dimensionless quantity. In general, the 
evaluation of Cj^^ is a bit cumbersome since there is 
no unique ground state and degenerate perturbation the- 
ory must be applied. When both A^i and N2 correspond 
to closed shells, then Cjt can be calculated straight- 
forwardly analytically [2l| , 

C^^^^^ = inaL J pr{r)p^Hr)df. (14) 

Here, p^^{r) is the density of a one-component non- 
interacting gas with Ni fermions of mass m^, normalized 
so that J ^ {r)dr = Ni. Alternatively, one can approxi- 
mate the pf^ by the Thomas- Fermi density profiles. This 
approximation should be quite accurate in the large N 
limit. 




5 10 15 20 

N 



FIG. 1: (Color Online) Cj^rj^jv^ coefficients divided by Eni as 
a function of A^. Circles correspond to L = ground state, 
squares to L = 1 ground state and triangles to L = 2 ground 
state. A solid line connects the odd-A values while a dashed 
line connects the even-A values. 

To obtain the Cfq^ for open-shell systems, we apply 
first-order degenerate perturbation theory. This calcula- 



TABLE I: Angular momentum L and coefficient Clf-^^pf^ for 
the ground state of equal-mass two-component Fermi gases in 
the weakly-attractive regime. Here, we consider A2 — Ni for 
even A and Ai = A2 + 1 for odd A. 
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tion additionally allows us to obtain the angular momen- 
tum quantum number L of the ground state. Figure [1] 
and Table [J present the results for A^ < 20 and k = 1. 
The coefficients increase monotonically with in- 

creasing A^ and show a slight odd-even staggering. In 
general, the coefficients Clf_^ for even N are compar- 
atively higher than those for odd N, implying a smaller 
energy for even A^ than for odd A^ and suggesting that, 
even in the perturbative regime, the odd-even oscillations 
are already present. We note that the coefficients 
for even A^ shown in Fig. [T] clearly reflect the shell-closure 
at A^ = 8. 

In the weakly-interacting molecular BEC regime, the 
two-component Fermi system should behave like a system 
that consists of Nd bosonic molecules and Nf = or 1 
fermions. In first order perturbation theory the ground 
state energy of such a system is given by 

p AT T^n 1^ I 3A^/ , Nd{Nd-l) [2 add 
E^NdE{l,l) + fu.^ + h^ ^--^ 

"■ho 

+n.NdN,^^il5) 

"ho 

Here, add and aad denote the dimer-dimer and atom- 
dimer scattering lengths, respectively. The oscilla- 
tor lengths ajjo'^'* and a^/^^^ for the dimer-dimer and 
atom-dimer systems, a^/ff^ — ^ h/ {2pdS^) and a'^y^f' — 
■y/ h/ {2padi^), are defined in terms of the reduced mass 
Pdd of the dimer-dimer system and the reduced mass pad 
of the atom-dimer system, respectively. 

The limiting behaviors of the BEC-BCS crossover 
curve can be used to guide the construction of the many- 
body nodal surface, which is a crucial ingredient for our 
FN-DMC calculations (see Sec. IhTB]) . In the weakly- 
interacting molecular BEC regime, even A'^ systems con- 
sist of A^/2 dimers. Each molecule is expected to be in 
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its rotational ground state, leading to a many-body wave 
function with total angular momentum L — 0. For odd 
iV systems, the extra fermion is expected to occupy the 
lowest s-wave orbital, leading, as in the even N case, to 
a many-body wave function with L = 0. Thus, the angu- 
lar momentum of even N systems is expected to be the 
same along the crossover while that of odd iV systems is 
expected to change (for N = 3, this has been pointed out 
recently by two independent groups [3l|, [111). This sym- 
metry change introduces a kink in the normalized energy 
curve A^^^-* for odd N and in the excitation gap A(iV) 
(see below) at the scattering length where the symmetry 
change or inversion occurs. 

In addition to the energy crossover curve, we calcu- 
late the excitation gap A(A^), which characterizes the 
odd-even oscillations of two-component Fermi systems, 
as a function of N. For homogeneous two-component 
Fermi systems with equal masses, the excitation gap 
A, which equals half the energy it takes to break a 
pair, is quite well understood. In the weakly interact- 
ing BCS regime, the excitation gap A becomes expo- 
nentially small [33I [3^ , indicating vanishingiy little pair- 
ing. In the deep BEC regime, on the other hand, the 
excitation gap approaches half the binding energy of the 
free-space dimer, indicating essentially complete pairing: 
By adding an extra particle to the odd N system, the 
energy of the total system changes by approximately 
the binding energy of the free-space dimer. In addi- 
tion to these limiting cases, the excitation gap of the 
equal-mass two-component Fermi system has been deter- 
mined throughout the crossover regime by the FN-DMC 
method [ssl IsgI [37|. For unequal-mass systems, in con- 
trast, the behavior of the gap is much less studied and 
understood [H [H, [H lU, El . 

To define the excitation gap A(7V) for trapped 
unequal- mass systems, we set iV = 2n -I- 1 and assume N 
to be odd. The unequal-mass system is characterized by 
two chemical potentials, the chemical potential /ii(iV) for 
species one and the chemical potential ^2 {N) for species 
two (see, e.g., Ref. HI), 

E{n + 1, n) = E{n, n) + ^i(2n + 1) + A(2n + 1) (16) 

and 

E{n, n + l) = E{n, n) + /i2(2n + 1) + A(2n + 1). (17) 

Here, A(2ri + 1) denotes the excitation gap. If A(2n -|- 
1) vanishes — as is the case for the normal system — , 
then Eqs. I|16p and p7)) reduce to the "usual" chemical 
potentials. Furthermore, Hi{N) and i^2{N) coincide for 
equal-mass systems. To determine /ii(iV), /i2(iV) and 
A(A'^), we need an additional relationship. In condensed 
matter physics, one typically considers the average of the 
two chemical potentials, 



[^ll{2n + 1) + fi2i2n + 1)] = 
\ [E{n + l,n+l) - E{n,n) 



Since the average chemical potential is defined in terms 
of the energy of the next smaller and the next larger 
balanced systems, it is independent of the odd-even os- 
cillations. Equations through (fTS]) can be solved for 
^i(2n -I- 1), ^2(2n + 1) and A(2n + 1), 



^il{2n + l) 



E{n + 1, n + 1) - E{n, n + 1) 
2 

E{n+ l,n) - E{n,n) 



, (19) 



E{n + l,n + l)-Ein + l,n) 
fi2[2n + 1) = 

E(n,n+ 1) - E{n,n) 



, (20) 



and 



A(2n-f 1) 



E{n + 1,71) + E{n,n + 1) 
' 2 
E{n,n) + E{n + l,n + l) 



(21) 



Note that the energies E{n + 1, n) and E{n, n + 1) are 
equal for equal masses. The excitation gap A(A^) and 
the chemical potentials fii{N) and fJ.2{N) depend on N, 
K, uj and fls. 

Ultimately, one of the goals is to relate the excitation 
gaps calculated for the trapped and the homogeneous 
systems. For equal masses and equal frequencies, the 
densities of the two trapped species overlap fully. Hence, 
one might expect that the excitation gaps of the homoge- 
neous and inhomogeneous systems can be related via the 
local density approximation (LDA), which predicts that 
A(A^) scales with N as N-^^^. Connecting the excitation 
gaps for the homogeneous and trapped systems in this 
way breaks down, however, if the extra particle sits near 
the edge of the gas cloud; this is the region that is poorly 
described by the LDA. Indeed, we present some evidence 
that the extra particle sits for > 11 near the cloud 
edge. For unequal masses, the connection between the 
two excitation gaps becomes even more challenging, be- 
cause one now has to first determine whether the trapped 
system exhibits phase separation or not [13, SI] . 

For the trapped system, the density mismatch can be 
quantified by comparing the density overlap OJvi JV2 ' 



(22) 



(18) 



of the unequal-mass system with that of the equal-mass 
system for a given scattering length Os. In Eq. (|22p . 
the one-body densities Pi{r) and the oscillator length 
ttho depend on n. In the non-interacting limit, the 
normalized density mismatch ATa/'^TVi reduces to 
(^NuN2/^Ni,N2- In t^^is case, O^^ ffJOjf^j^^ equals one 
for all K if TVi = A^a = 1 (see Table II of Ref. (U). For 
larger A^, however, OJ^^ n2/^Ni N2 decreases from 1 to 
a finite value that is smaller than one as k varies from 
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one to infinity. In particular, the Thomas Fermi approxi- 
mation predicts 0%^^j^^/0]^^^^^ 3157r/1024V2 w 0.683 
for large non- interacting systems (A^i — N2) with large 
AC. For the small unequal-mass systems considered in 
Sec. llVi we find that the density mismatch for finite 
is smaller than that for 0= = 0. 



C. Hyperspherical formulation at unitarity 

The two-component Fermi gas at unitarity is charac- 
terized by a diverging scattering length, i.e., l/a^ — 0. In 
this regime, the underlying two-body potential, for suffi- 
ciently small i?o, has no characteristic length scale, thus 
leaving only the size of the system itself. This elimina- 
tion of the two-body length scale is the key to obtain- 
ing a number of analytical results; a particularly appeal- 
ing framework for deriving these results employs hyper- 
spherical coordinates. The hyperspherical formulation 
has been primarily develo ped in the context of few-body 
systems [Isl. lielliTl. lislligl. l50l| . More recently, some prop- 
erties of Bose and Fermi gases with essentially arbitrary 
number of atoms have been explained successfully within 
this formulation [5l|, [13, US] ■ The ability to treat both 
small and large systems on equal footing makes the hy- 
perspherical formulation particularly suited for studying 
the transition from few- to many-body systems. 

We define the hyperspherical coordinates by first sep- 
arating off the center-of-niass vector Rcm, and by then 
dividing the remaining 3N — 3 coordinates into the hy- 
perradius R and 37V — 4 hyperangles, collectively denoted 
by Q. The hyperradius R is defined by 
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(23) 



and can be viewed as a coordinate that measures the 
overall size of the system. Here, M denotes the total 
mass of the system, M = miNi + m2N2, and an 
arbitrary mass scaling factor. Usually, the value of /i^v is 
chosen so that the hyperradial potential curves Vs^{R), 
defined below, approach physically motivated asymptotic 
values as i? — > 00. 

In the adiabatic approximation [4^ . the relative wave 
fimction ^''''^'(i?, il) reduces to 

r'^'(i?,r!) = R-^^^-^'^/^F,niR)'^AR;n). (24) 

The antisymmetric Pauli correlations are built into the 
channel functions ^^(i?;^!) at the outset. In addition, 
the $i/(i?; fl) account for a significant fraction of the two- 
body correlations of the system. Within the hyperspher- 
ical approximation, the description of the many-body 
system reduces to solving a one-dimensional Schrodinger 
equation in the hyperradial coordinate R, 



1 



2^7V dR^ 



The effective hyperradial potential Vs^ (R) includes part 
of the kinetic energy and a contribution due to the short- 
range two-body interactions. 

Assuming zero-range interactions, the adiabatic ap- 
proximation becomes exact for a subclass of univer- 
sal states of the unitary two-component Fermi gas [H^ . 
For these states, the channel functions obey specific 
boundary conditions imposed by the zero-range pseu- 
dopotential and become independent of R. Furthermore, 
the functional form of the hyperradial potentials 14,^ {R) 
can be derived analytically |5J, [sll , 



VsAR) 



h^s^{si, + 1) 



2^Ari?2 

The eigen energies of Eq. ([^5]) are then given by 



rprel 



2n 



— 1 

2/ ' 



(26) 



(27) 



where n is a non- negative integer, and the hyperradial 
wave functions F^n{R) (not normalized) by 



2C^ 



(28) 



EltF,^{R). (25) 



where C denotes the oscillator length associated with /xjv, 
C — y/K/JJiMUj) , and Ln"'^'^^'^^ the Laguerre polynomial. 
The total energy Ei,n is obtained from by adding 
the center of mass energy. The spacing between states 
labeled by the same v is 2hLij and is thus independent of 
Sy. This implies that knowledge of the lowest eigenenergy 
in each hyperradial potential curve determines the 
entire energy spectrum. This property of the spectrum 
has also been shown using the scale invariance properties 
of unitary systems [56j . Transitions between vibrational 
levels that lie within a given hyperradial potential curve 
Vs„ (R) can be driven by an excitation operator that de- 
pends on R only. Such a driving field results in a ladder 
of excitation frequencies of the form 2khijj, where k de- 
notes an integer. On the other hand, transitions between 
states living in different hyperradial potential curves (la- 
beled by v and v') require the driving field to depend 
on Q, or stated more generally, the excitation operator 
must not commute with the fixed-hyperradius Hamilto- 
nian. The corresponding excitation frequencies are, in 
general, non-integer multiples of 2hLU and depend on the 
difference between Si, and Si^/. Thus, knowledge of the 
entire excitation spectrum requires determining all s^. 
Moreover, the coefficients of the three-body system 
play a role in determining the three-body recombination 
rate for large and negative Os (57| . and the lifetime of 
weakly bound dimers for large and positive as 57]. Sim- 
ilarly, one may expect that the of larger systems play 
a role in determining the corresponding quantities for 
larger systems. Section HV CI presents evidence of the 2huj 
energy spacing and determines the coefficients for the 
four-particle system for various mass ratios. 
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Equation (|23p defines the hyperradius R without the 
CM motion. Alternatively, we can define a hyperradius 
R', 



MR'^ = UnR^+MR, 



(29) 



which includes the CM motion and represents the rms 
radius of the system. In the adiabatic approxima- 
tion, the total wave function ^'(i?',f2') can be writ- 
ten in terms of the new hyperradius R' as ^(i?', ri') — 
i?'-(3^-i)/2F^„(i?')$(i?';rj'), where fl' collectively de- 
notes the 3A^ — 1 hyperangles. Equations (pS)) and 
remain valid if R, fiN and F^n are replaced by i?', M and 
F^n, respectively. The eigen values of the hyperradial 
Schrodinger equation equal the eigenenergies E^in of the 
total system. Defining x = R' /R'j^j and e^n = E^n/Eiqi, 
the hyperradial Schrodinger equation can be rewritten as 



1 



2fleff dx^ 2fieffX^ 



= ^unFyn{x), (30) 



where /^e// = Ej^j/{hLj) . Above, R'j^j denotes the rms 
radius of the non-interacting system; it can be, using the 
virial theorem (s^ . [ssj , expressed in terms of the energy 
Eni of the non-interacting two-component Fermi gas. 



h I Eni 



The dimensionless coefhcients C 



c 



N 



so(so + l) so{so + l)h'^uj^ 



E^ 



(31) 



(32) 



characterize the ground state of the system at unitarity. 
The scaled hyperradius x and the scaled energies e^n re- 
main finite in the large N limit and are thus particularly 
well suited to discuss the large N limit (see Sec. I IV Bp . 
For small systems, in contrast, some properties of the sys- 
tem can be highlighted more naturally using the unsealed 
hyperradius R or R' . 

The coefficients describe both the trapped and free 
systems, and can be related to the universal parameter 
^ of the homogeneous system The hyperspherical 

framework thus connects few- and many-body quantities 
and allows one to bridge the gap between atomic and 
condensed matter physics. 



III. NUMERICAL TECHNIQUES 
A. Correlated Gaussian approach 

The CG method has proven capable of providing an 
accurate description of trapped few-body systems with 
short-range interactions [lo, [SI, ill. The CG method 



expands the many-body wave function ^' in terms of a 
set of basis functions ^^^^.^-j, 

*(ri,-- - ,fN)^ C{d,,}*{d.,}(ri,-- - ,fAr), (33) 

{dij} 

where the C^dij} denote expansion coefhcients and the 
{dij} a set of widths. Each basis function has the form: 



{ V'o(-RcJ\/) exp 




(34) 

Here, -00 is the ground state wavefunction associated with 
the center-of-mass vector Rcm, and the operator S en- 
sures that the basis functions have the proper symmetry 
under exchange of two fermions of the same species. Due 
to the simplicity of the basis functions, the elements of 
the Hamiltonian and overlap matrices can be calculated 
analytically [g^, [6l|- Since the basis functions depend 
only on the center of mass vector and the interparticle 
distances, i.e., Gaussians centered around = 0, the re- 
sulting eigenenergies correspond to eigenstates with zero 
relative angular momentum Lrei and zero total angular 
momentum L; throughout this work, we do not consider 
center-of-mass excitations so that Lrd = L for all sys- 
tems investigated. To determine the eigenenergies of 
states of the iV-atom system with non- vanishing Lrei, we 
add a spectator atom and solve the Schrodinger equation 
for the (iV-t- l)-atom system. The extra particle does not 
interact with the rest of the system but can have non- 
vanishing angular momentum. This trick allows us to 
describe non-zero angular momentum states of the N- 
atom system. We find, e.g., that the ground state of the 
equal-mass three- and five-particle systems at unitarity 
has Lrei = 1- 

To illustrate how the energies calculated by the CG 
method converge with respect to the size of the basis set, 
we consider the three-body system with L = at uni- 
tarity. We define Ed as the eigenenergies obtained for 
an optimized basis set of size D. The optimization of 
the basis functions for a given size D is performed us- 
ing the basic ideas of the stochastic variational approach 
6l|. The size of the basis set is then increased and the 
new basis functions are optimized. Figure [5] shows an ex- 
ample of the convergence of the lowest few eigenenergies 
for Rq = O.Olaho as a function of D. The largest D con- 
sidered in this study is 700, and the energies have been 
tested and are approximately converged for this D value. 
Thus, Fig. [5] shows the normalized difference between 
Ejoo and Ed for the lowest few eigenenergies. Figured] 
shows that the basis set can be improved systematically. 

For larger number of particles, the size of the basis 
set needs to be increased. For N — 5 and 6, the size of 
the basis set is increased up to approximately D = 10^. 
The = 6 energies reported in Ref. [l^l, e.g., are cal- 
culated for D = 1.6 X 10^. Here, we analyze the con- 
vergence of these energies as a function of 1/D. Since 
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100 200 300 400 500 600 

D 

FIG. 2: Convergence of the energetically lowest-lying energies 
as a function of the size D of the basis set for N = 3 {L = 0) 
at unitarity. The range _Ro is fixed at O.Olaho- Solid lines 
connect the CG energies (filled circles) of a given state for 
ease of viewing. 



the energies behave approximately hnearly as a func- 
tion 1/D, we can extrapolate straightforwardly to the 
limit _D — > oo. The extrapolated energies for v = 
are Eqo = 8.48ftw, Eqi = lO.bOhoj and £^02 = 12.50huj. 
Eqo and Eqi agree with those reported in Ref. |23| for 
D ~ 1.6 X while £^02 is only Q.02huj lower than the 
previously reported value. For v = 1 and 2, the extrap- 
olated energies are £10 = lOASfuv and £20 = 10-99hiu; 
these energies are lower by O.Olhui than those reported 
in Ref. [20|. While the extrapolated energies are most 
likely closer to the exact eigenenergies than the energies 
calculated for D = 1.6 x lO**, we note that the extrap- 
olated energies are no longer variational, i.e., they no 
longer provide upper bounds to the exact eigenenergies. 
Our analysis of the v < 2 excited energies shows that 
the extrapolated energies follow the expected 2huj spac- 
ing more closely than those calculated for the largest D 
considered, suggesting that the extrapolation procedure 
is indeed justified. 

In general, the convergence of the energies with re- 
spect to the basis set depends on the scattering length 
Gs and the number of states considered. Usually, an accu- 
rate determination of the spectrum at unitarity requires 
a larger basis than the determination of the spectrum 
on both the weakly- interacting BEG and BCS sides. For 
equal-mass systems, a converged basis at unitarity usu- 
ally describes the spectrum in the entire crossover re- 
gion accurately. Of the equal-mass systems treated, the 
= 5 (L = 1) calculations have been the hardest to 
converge. For L = 1 states, we can estimate the un- 
certainty of the calculations by monitoring the energy of 
the spare non-interacting particle, which is known ana- 
lytically. For example, for the N — 5 equal-mass calcu- 
lations presented in Ref. [lO], the energies of the spare 
non-interacting particle deviate from the exact solution 
by approximately O.Olhuj, which is less than 1%. We find 



that systems with large k are typically harder to converge 
than the corresponding equal-mass systems. 

To analyze the effects of finite range interactions 
we study the eigenenergies of the three-particle system 
at unitarity (L = 0) as a function of the range Rq- 
Figures ^a) through (c) show the energies for the low- 
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FIG. 3: (Color Online) Three-body energies E^o at unitarity 
for L = [(a) v = Q, (b) v — 1 and (c) 1/ = 2] as a, function of 
the range Ro- Symbols show the CG energies and solid lines 
the hnear extrapolation to the _Ro = limit. 

est state in the hyperradial potential curve Vs^ (i?) with 
J' = 0, 1 and 2. The energies show a linear dependence 
on i?o, and can thus be extrapolated straightforwardly 
to the zero range limit. Neglecting the basis set er- 
ror, which is estimated to be smaller than the uncer- 
tainty of the extrapolation, we find Eqo — 4.66622(l)/ia;, 
£10 = 7.62738(2)ftw, and £20 = 9.61466(4)ftu;. Our 
three-body energies compare favorably with those cal- 
culated using the Si, coefficients, v — and 1, deter- 
mined by Ref. [57] in Eq. ([27]), ^q q = 4.6662220fiw and 
Eio = 7.6273521/iw. Section flV Al reports three-particle 
energies for equal masses for Rq — O.Olaho, which — ac- 
cording to Fig. [3] — agree to better than 0.02ftti; with 
those calculated in the zero-range limit. We addition- 
ally performed systematic studies of the dependence of 
the energies on the range Rq for the three-body system 
with equal and unequal masses in the weakly-interacting 
molecular BEC regime, where two-body bound states 
form (see Sees. Ill Bl and IIV A|) . and for the four-body 
system. For the five- and six-body calculations, it is pro- 
hibitively expensive to perform calculations for different 
Rq. For these systems, we estimate the finite range ef- 
fects based on our findings for the N = 3 and 4 systems. 

In addition to the energies. Sec. IIVDI reports struc- 
tural properties calculated by the CG approach. The 
one-body density and the pair-distribution functions are 
extracted from the total wave function 5* calculated by 
the GG approach by integrating 'i>^ over the relevant Ja- 
cobi coordinates. 
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B. Fixed- node diffusion Monte Carlo approach 

For larger systems, the CG approach in our current 
implementation becomes prohibitively expensive and we 
instead determine first-principles solutions of the time- 
independent Schrodinger equation using Monte Carlo 
techniques. 

In this study, we use the FN-DMC method [H ill, 
a variant of the diffusion Monte Carlo (DMC) method, 
to determine solutions for up to = 30 fermions. The 
DMC method, which interprets the system's wave func- 
tion as a density, allows for the accurate determination 
of the energy of nodeless ground states but is not suited 
to determine the energy of excited states of bosonic sys- 
tems or of fermionic systems. To treat systems whose 
eigenfunctions have nodes, the DMC algorithm has to be 
modified slightly. Here, we adopt the FN-DMC method, 
which obtains a solution of the Schrodinger equation that 
has the same symmetry as a so-called guiding function 
ipT- The FN-DMC method provides, to within statistical 
uncertainties, an upper bound to the exact eigen energy 
of the many-boson or many-fermion system, i.e., to the 
lowest-lying state with the same symmetry as ipT- 

If the nodal surface of ipx coincides with that of the 
exact eigenf unction, then the FN-DMC method results in 
the exact eigen energy of the system. In general, however, 
the nodal surface of the exact eigenfunction is not known 
and the FN-DMC results depend crucially on the quality 
of the nodal surface of ^t- In this work, we consider 
three different parametrizations of the nodal surface of 
two-component Fermi systems. 

The guiding function ipTi reads 

i=l i' = l 

The function F^^^^ determines the nodal structure ofipri 
and is, for even N and A^i = N2, constructed by anti- 
symmetrizing a product of pair functions / (6^ . 

Fnode^Afirii')J{r22')r--JirmN2)), (36) 

where A is the antisymmetrization operator. The pair 
function / is given by the free-space two-body solu- 
tion [g^] : / coincides with the free-space two-body bound 
state solution for positive scattering length a^, and with 
the free-space scattering solution, calculated at the scat- 
tering energy Erei, for negative as- For A'^ = 6, we 
treat Erei as a variational parameter and find a reduc- 
tion of the energy of 1 or 2% for a finite Erei compared 
to Erei = 0. For larger N, we simply use Erei = 0. For 
odd iV, we add a single particle orbital 0„; in Eq. ([55]) so 



that becomes, ioi Ni = N2 + I [Hiil, 

■ - ■ 

Aifinv),--- : .firNi-iM2), (t^niirNjal^^)) = 
finv) ■■■ /(riArJ <pniiri/a'f^^) 
/(?-2i') ••• f{r2N2) 0n;(^2/aio) 

det 

f{rN^l') ■■■ f{rN,N2) (t>nl{rN^/a\^^) 



where aj^^ = ^h/(miUj). We consider a number of differ- 
ent single particle orbitals <j)ni , and determine the optimal 
nl values by performing a series of FN-DMC calculations. 
For the lowest n and Z, the orbitals read (j^wi'P / o-^hl) — 1; 

0oi(r7ait^) = z/a«, </>2o{r7ait^) = 1 - 2(r/a«)V3 and 
0O2(r7ait^)=3(z/a«)2-(r/aW)2. 

In Eq. the $i (i = 1 and 2) denote Gaussian 

single particle orbitals that depend on a width parame- 
ter 6„ $,(rO = exp(-rV(26f)). If h = ^h/{m,Lo), 
coincides with the ground state orbital of the harmonic 
oscillator. The parameters bi and &2 are optimized vari- 
ationally. For even N {Ni = N2) and equal masses, we 
require bi = 62. At unitarity, e.g., we find that the bi 
are smaller than the , reflecting the attractive nature 
of the interspecies interaction potential. If hi = , the 
product ^i{r)4'ni(f/ afl) equals the harmonic oscillator 

wave function 4>^rda \r / a(l) . 

In Eq. ([55)l . the pair functions gn, 322 and 512 are 
introduced to improve the variational energy and to ad- 
ditionally ensure that the structural properties calculated 
at the VMC and FN-DMC levels agree at least qualita- 
tively. The pair functions 51 1 and (722 allow for the ef- 
fective repulsion between equal fermions to be accounted 
for, 

gii{r) = exp{~pir~'") (38) 

for i — I and 2. The parameters pi, p2, qi and 52 are 
optimized variationally. For even A^ and equal masses, 
we require pi — p2 and qi = q2- The pair function gi2 is 
parametrized in terms of the three variational parameters 
t, P12 and gi2, 

gi2{r) - 1 + texp(-pi2r-«i=). (39) 

The parameters t, pi2 and qi2 are optimized under the 
constrained that gi2 > 0. 

The guiding function is expected to provide a good 
description of the system in the weakly-interacting molec- 
ular BEG regime, where we expect bound dimer pairs to 
form. Section IIVBI shows that this wave function also 
provides a good description of the unitary gas for suf- 
ficiently large N. This is in agreement with FN-DMC 
studies for the homogeneous system (6^ . Since each pair 
function / has vanishing relative orbital angular momen- 
tum, the total angular momentum L of ^pTi is for even 
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N and Ni = N2- For odd N, L of V'Ti is determined by 
the angular momentum of 0„;, i.e., L — I. 

In addition to ipxi, we consider the guiding function 

Ni N2 

i^T2 = n*i(^*) ^ n ^2(r»') X 

i=l i' = l 

<ode(^i,-- - ,r^Jx n /(r,,0- (40) 

The nodal surface of ipT2 is determined by '^nodc which is 
defined so that the product Ylf^i x Yl^ti ^2{ri') x 

^node coincides for bi = a^j^^ with the wave function of 
N trapped non-interacting fermions. Thus, the nodal 
surface of '(/'T2 coincides with that of the corresponding 
non-interacting system. The pair function / coincides 
with the pair function / introduced above for r < Rm, 
where Rm is a matching point determined variationally. 
For r > Rm, / is given by ci -I- C2exp(— ar). The pa- 
rameters ci and C2 are determined by the condition that 
/ and its derivative be continuous at r = Rm while a is 
optimized variationally. 

The guiding function '^t2 is expected to provide a good 
description of the system in the weakly-interacting BCS 
regime. In this regime, we construct the guiding function 
so that its angular momentum agrees with that predicted 
analytically (see Table [l|. Section FlV BI shows that the 
guiding function '^t2 also provides a good description of 
small fermionic systems at unitarity. 

Finally, the guiding function 'ijjxa is constructed follow- 
ing Eqs. (3) and (4) of Ref. [H. We find that Vts gives 
the lowest energy for TV = 11. 

Expectation values (A) of operators A that do not 
commute with the Hamiltonian cannot be calculated as 
straightforwardly by the FN-DMC method as the energy. 
Here, we use the mixed estimator {A)mixed [63l. Igt}. 

{A)rmxed = 2{A)dMC - {A)vMC- (41) 

In Eq. ([1T|) . {A)vMC denotes the expectation value calcu- 
lated by the VMC method and {A) d mc that calculated 
by the FN-DMC method. We note that some algorithms 
for the calculation of pure estimators exist [sl, [6§| but 
we do not use them in this work. 

In some cases, we optimize the variational parameters, 
collectively denoted by p, by not only minimizing the en- 
ergy expectation value but by additionally ensuring that 
ipT captures selected structural properties correctly. To 
this end, we compare the structural properties calculated 
by the VMC method for a given po with those obtained 
by the FN-DMC method, which uses ipripo) as a guid- 
ing function, and then choose a new parameter set pi 
so that the VMC structural properties calculated using 
Pi agree better with the FN-DMC structural properties 
calculated using po- This procedure is repeated till the 
VMC and FN-DMC structural properties and energy ex- 
pectation values agree sufhciently well. For equal-mass 



TABLE II: Dimer-dimer scattering length add and dimer- 
dimer effective range Vdd obtained using (a) the CG spectrum 
and (b) the FN-DMC energies. The reported uncertainties 
reflect the uncertainties due to the fitting procedure; the po- 
tential limitations of the FN-DMC method to accurately de- 
scribe the energetically lowest-lying gas-like state, e.g., are 
not included here (see Sec. IIIB of Rof. 21]). 



K 


add /as (a) 


add/as (b) 


rdd/as (a) 


Tdd/as (b) 


1 


0.608(2) 


0.64(1) 


0.13(2) 


0.12(4) 


4 


0.77(1) 


0.79(1) 


0.15(1) 


0.23(1) 


8 


0.96(1) 


0.98(1) 


0.28(1) 


0.38(2) 


12 


1.10(1) 


1.08(2) 


0.39(2) 


0.55(2) 


16 


1.20(1) 


1.21(3) 


0.55(2) 


0.60(5) 


20 


1.27(2) 


1.26(5) 


0.68(2) 


0.74(5) 



systems with N < 20, our VMC energies are at most 
15% higher than the corresponding FN-DMC energies. 
The optimization strategy em ploy ed here is similar in 
spirit to that discussed in Ref. [70] for the homogeneous 
system. 



IV. RESULTS 

A. Ground state energy in the crossover regime 

This section discusses the behavior of the crossover 
curve and the excitation gap for = 3 for different mass 
ratios k. This odd N study complements our earlier re- 
sults for even N [2l[. Our analysis for = 4 showed 

that the crossover curve A^^-* is independent of the de- 
tails of the two-body potential and allowed us to extract 
the dimer-dimer scattering length add and the dimer- 
dimer effective range rdd as a function of n. The a^d and 
Vdd results from Ref. '2l[ are summarized in Table [TTl 
Furthermore, for larger even N systems, we determined 
the validity regimes of the analytically calculated limit- 
ing behaviors in the weakly-interacting molecular BEC 
and atomic BCS regimes. Our even N study resulted 
in a deeper understanding of some of the peculiarities of 
trapped systems and emphasized similarities and differ- 
ences between the trapped and homogeneous systems. 

The behavior of odd N systems is rich and, in many 
cases, qualitatively different from that of even N sys- 
tems. One characteristic of odd N systems is the pos- 
sible change of the angular momentum of the ground 
state as the scattering length is tuned through the BEC- 
BCS crossover region (see Sec. IIIBI and Refs. [3l|, Is^). 
Figure [3] shows the three-particle energy E, with the en- 
ergy E{1, 1) + 3hijj/2 subtracted, for L = (solid lines) 
and L — 1 (dashed lines). The upper panel shows results 
for K, — I, and the two lower panels for k = 4 Jpan- 
els (b) and (c) consider the three-particle system with a 
spare heavy and a spare light particle, respectively]. The 
ground state has L = 1 for aho/a-s —oo and L = 
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FIG. 4: (Color Online) Normalized energy {E - E{1, 1) - 
?>t\w /2) /huj for = 3 as a function of aho/a-s calculated by 
the CG approach (lines). E denotes the three-body energy 
for L = (solid lines) and for L = 1 (dashed lines): (a) 
equal-mass atoms [k = 1, i5 = E{2,1) = E{1,2)\, (b) two 
heavy atoms and one light atom [k = 4, E — E{2,1)], and 
(c) two light atoms and one heavy atom [k = 4, E = E{1, 2)]. 
The normalized energy crossover curve Ag''' , Eq. ((5} , coincides 
with the dashed and solid lines, respectively, depending on 
whether the three-particle ground state has L = 1 or 0. In 
the GG calculations, the range _Ro of the two-body potential 
is fixed at O.Olaho- For comparison, crosses and circles show 
selected FN-DMC energies for L = and L = 1, respectively. 



for aho/o-s — ^ oo, independent of k and independent of 
whether the spare particle is heavy or light. For equal 
masses, the change of symmetry occurs at « aho- For 
K = 4, in contrast, it occurs at ~ O.Saho if the extra 
particle is a heavy atom [panel (b)] and at Ug ~ 3aho if 
the extra particle is a light atom [panel (c)]. The dashed 
and solid lines shown in Fig. 3] coincide with the normal- 
ized crossover curve A^^' Eq. ([5]), in the region where 
the ground state of the three-particle system has L = 1 
and 0, respectively. The normalized crossover curve Ag'''' 
changes from 1 in the weakly-interacting molecular BEC 
regime to in the weakly-interacting BCS regime. 

We find that the normalized L = 1 energy curve for 
two heavy atoms and one light atom [Fig. |D^b)] depends 
notably on the range of the underlying two-body poten- 
tial if the scattering length Os is positive. For example, 
the normalized energy curve changes by as much as 20% 
if the range Rq of the two-body potential changes from 
O.Olaho to 0.02a/io. This comparatively large dependence 
on Rq indicates that the properties of the system with two 
heavy atoms and one light atom are not fully determined 
by the s-wave scattering length for the ranges consid- 
ered. In the -Ro ^ limit, the k = 4 system is expected 
to behave universal [H, [23| . We speculate that the com- 
paratively strong dependence of the normalized energy 
curve on the range for > is related to the fact that 
the three-particle system supports, for sufficiently large 



AC, bound states with negative energy. 

For comparison, circles and crosses in Fig. [3] show se- 
lected three-particle energies calculated by the FN-DMC 
method for L = and L = 1, respectively. The good 
agreement with the CG results (lines) indicates that the 
FN-DMC method can be used to accurately describe dif- 
ferent symmetry states. 

Our CG energies for equal-mass systems interacting 
through short-range potentials presented in Fig. CTa) can 
be compared with those of Kestner and Duan [31[ ob- 
tained for zero-range interactions. Our L — 1 energy 
curve agrees with that of Kestner and Duan for all scat- 
tering lengths as considered. The L — energy curve, 
however, only agrees for as < 0. For as > 0, our results 
are noticeably lower than those of Kestner and Duan. 
As shown below, our > results for L = predict 
the correct atom-dimer scattering length suggesting that 
our energies should be very close to those for R^ — 
and that the disagreement is not due to finite-range ef- 
fects. We speculate that the results of Kestner and Duan 
might not be fully converged for > although other 
possibilities cannot be excluded. 

Figures \5la) and (b) present the BCS and BEC lim- 
iting behaviors for an equal mass system with iV = 3. 
The perturbative expression, Eq. (fT5|) . on the BEC side 
is expected to be applicable if i?o ^ *C o-ho', thus, 
we choose a small Rq, i.e., Rq — O.OObaho, in the CG 
calculations. The energy is in this region determined by 
the atom-dimer scattering length Uad [see Eq. (fTS)) ]. The 
CG energies change linearly with a^., showing that aad 
is proportional to a^., i.e., aad = CadOis- A simple lin- 
ear fit to the CG results predicts Cad ~ 1-21, in good 
agreement with previous studies [tH . [t^ . which found 
Qad ~ l-2as- A solid line in Fig.lHlIa) shows the resulting 
linear expression. A more sophisticated analysis accounts 
for the energy-dependence of aad [zl, [ill , which results 
in a more reliable determination of Cnd and also a deter- 
mination of the effective range Tad [2ll |. Considering the 
three lowest energy levels on the BEC side [2l[, we ob- 
tain Ca d ~ 1.18(1) and Vad ~ 0.08(1)0^. It was suggested 
earlier |24| that the atom-dimer system is characterized 
by a soft-core repulsion with range of the order of a^; 
our calculations support this general picture but predict 
a range about ten times smaller than a^. On the BCS 
side, the first order correction varies also linearly with as- 
Circles in Fig. [Hb) show the CG results while the solid 
line shows the prediction from Eq. p3p . Good agreement 
is observed in both limiting behaviors. 

Our CG energies for N = 2, 3 and 4 can be readily 
combined to determine the excitation gap A (3), Eq. (|2ip . 
Figure [6] shows the excitation gap A(3) as a function of 
o-ho/o-s for two different mass ratios, i.e., k = 1 and 4. 
In the weakly-interacting molecular BEC regime, the ex- 
citation gap approaches 3?ia;/2 — E(l, l)/2 (circles), in- 
dependent of the mass ratio. In the weakly-interacting 
BCS regime, however, the excitation gap depends on the 
mass ratio (see inset of Fig. 15]). For equal masses, A(3) 
is very well described by the perturbative expression for 
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FIG. 5: (Color Online) Limiting behavior of the ground state 
energy for A'^ = 3 equal mass fermions. (a) Energy correction 
AE = E{2, 1) - E{1, 1) - 3hLj/2 on the BEC side. Circles 
show the CG results while the solid line shows the first order 
correction for Uad ~ 1.2as. (b) Energy E{2,1) on the BCS 
side. Circles show the CG results while the solid line shows 
the first order correction on the BCS side. 
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FIG. 6: (Color Onhne) Excitation gap A(iV) for iV = 3 as a 
function of ahojas calculated by the CG approach for k = \ 
(solid line) and «: = 4 (dashed line). Circles present the BEC 
limiting behavior ihiu/2 — E{1, l)/2 which is independent of 
K. The inset shows a blow-up of the region where A(3) is 
smallest; in this region, the dependence of A(3) on k is most 
pronounced. The dash-dotted line shows the limiting behav- 
ior for K = 1 obtained by approximating the E{N) in Eq. (|2H) 
by their perturbative values, Eq. (|13p. 



Qs < —Q.5aho (dash-dotted line in the inset). Figure [H] 
shows that A(3) is smaller for k = 4 than for k — 1. 
Intuitively, this might be expected since the radial den- 
sities of the two species do not fully overlap for unequal 
masses (recall, we consider the case where species one and 
two experience the same trapping frequency). Thus, the 
pairing mechanism is expected to be less efficient in the 
unequal-mass system, especially on the BCS side, than 
in the equal-mass system. The next section discusses the 



TABLE III: CG and FN-DMC energies E at unitarity for 
small equal-mass systems with angular momentum L = 
and 1. The CG energies are calculated for the Gaussian in- 
teraction potential with Rq — COla^o for A'' = 3 and 4, and 
with Ro = 0.05aho for iV = 5 and 6. The FN-DMC energies 
are calculated for the square well interaction potential with 
Ro = O.Olafto. The guiding functions ipri and tpT2, Eqs. psp 
and (|40p . are used to obtain the energies of states with L — 
and 1, respectively. 



iV 


L 


E/{hoj) (CG) 


E/{riLu) (FN-DMC) 


3 





4.682 


4.67(3) 


3 


1 


4.275 


4.281(4) 


4 





5.028 


5.051(9) 


5 





8.03 


8.10(3) 


5 


1 


7.53 


7.61(1) 


6 





8.48 


8.64(3) 


7 







11.85(5) 


7 


1 




11.36(2) 


8 







12.58(3) 


9 







15.84(6) 


9 


1 




15.69(1) 



behavior of the excitation gap at unitarity in more detail. 

B. Ground state energy at unitarity 

This section explores the odd-even behavior of two- 
component Fermi gases at unitarity. In particular, we 
present the excitation gap for equal-mass systems with 
up to iV = 30 fermions and interpret the behaviors of 
these systems within the hyperspherical framework. We 
also discuss the excitation gap for small unequal-mass 
systems. 

Table [m] summarizes selected CG and FN-DMC en- 
ergies for small equal-mass systems at unitarity. Some 
of these energies were already reported in Refs. [l^, |2l| , 
and we include them in Table IIIII for comparative pur- 
poses. A comparison of the CG and FN-DMC energies 
for < 6 shows that the FN-DMC energies agree to 
within 2% with the CG energies for both L = and 
1 states. This agreement suggests that the nodal sur- 
face used in the FN-DMC calculations is quite accurate. 
Thus, Table [ml shows that the FN-DMC method allows 
not only for an accurate description of the ground state 
but also of excited states. For = 9, the energy of the 
L — 1 state is by only about 0.15?ia; smaller than that 
of the L — Q state. The ground state energies for larger 
N are reported in Table II of Ref. [13] • For both even 
and odd N (N > 9), we find that the angular momen- 
tum of the lowest energy state at unitarity is zero. Our 
FN-DMC energies thus suggest that the total angular 
momentum of the lowest energy states at unitarity has 
L = 1 for small odd A^ systems and L = for larger odd 
N systems. We note that this conclusion depends cru- 
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cially on the construction of the nodal surface entering 
the FN-DMC calculations. For N — 19, e.g., the energies 
at unitarity for L = 2 and 1 are less than O.Sfvjj higher 
than the i = energy; thus, the definite determination 
of the ordering of the states at unitarity with different 
angular momenta remains a challenge for odd-A^ systems 
with > 9. 

For homogeneous systems, the ground state energy per 
particle at unitarity £'„ is related to the energy per par- 
ticle Epc of the non-interacting system by a universal 
proportionality constant ^, = ^Epc [3l,[33,[63|- Ap- 
plying this result to the trapped unitary system with even 
N through the LDA, the ground state energy Eoo{N) of 
the trapped system becomes directly proportional to the 
energy E^i of the non-interacting trapped system [2lj |. 

EooiN) = V^Eni. (42) 

An analysis of our FN-DMC energies for = 2 — 30 
suggests that the trapped unitary system shows little 
shell structure. This motivates us to "smooth" the non- 
interacting energies, i.e., we approximate E^j by the ex- 
tended Thomas- Fermi expression (TSj . 

Eni,etf = n^i(3iV)4A-' (^1 + i(37V)-2/3^ . (43) 

To determine the proportionality constant ^, we fit our 
even N energies for TV = 2 — 30 to the expression 
V^Eni,etf- We find £,tr = 0.467, and denote the re- 
sulting energies by Efn . This value is in very good agree- 
ment with our previous result, S,tr — 0.465, obtained by 
including the energies for = 2 — 20 only [2l|. Cir- 
cles in Fig. [7] show the residual energy Eoo{N) — Efu for 
both even and odd A^. For even A^, the energy difference 
Eoo{N) — Efit is at most Q.lbhuj (except for A^ = 30, 
for which the error bar is large). This suggests that the 
energies of the trapped unitary system are indeed quite 
well described by \/itrENi.ETF] in other words, our en- 
ergies show little residual shell structure. As expected, 
the odd A^ energies are not even quantitatively described 
correctly by Eoo{N) — Efit. Instead, Fig.[7]shows that the 
residual energy Eoo{N) ~ Efu for odd A^ (circles) agrees 
quite well with the excitation gap A{N) (squares). For 
comparison, triangles in Fig. [7| show the excitation gap 
calculated using DFT [tI] . The good agreement between 
the DFT and FN-DMC results is encouraging. 

The ground-state energies Eoo{N) determine the co- 
efficients So [see Eq. (P7|) ] of the hyperradial potential 
Vso{R) [see Eq. (gB])]. Figures[5Ka) and (b) show the low- 
est hyperradial potential curves V{R) [V{R) ~ Vs„{R) + 
VtrapiR), where VtrapiR) = \^J,N^^'^R^ and = m] for 
A^ = 3 — 20 in the non-interacting limit and at unitarity, 
respectively. The small R behavior of V{R) is domi- 
nated by (R) while the large R behavior of V{R) is 
dominated by Vtrap{R)- Comparison of Figs. [DJa) and 
(b) shows that the attractive interactions lead to a low- 
ering of the potential curves at unitarity compared to 
those of the non-interacting system. Furthermore, the 
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FIG. 7: (Color Online) Excitation gap A(A'^) (squares) and 
residual energy Eoo{N) — Efu (circles) for equal-mass Fermi 
systems at unitarity as a function of A'' calculated from the 
FN-DMC energies. Triangles show A(A'^) calculated using 
density functional theory [761 ]. 
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FIG. 8: Hyperradial potential curves V{R) for equal-mass 
two-component Fermi systems with (a) vanishing interactions 
and (b) infinitely strong interactions as a function of R. The 
hyperradial potential curves naturally appear ordered as A'' 
increases: Solid lines correspond, from bottom to top, to = 
4 — 20 (A'^ even), while dashed lines correspond, from bottom 
to top, to TV = 3 - 19 (Af odd). 



V{R) at unitarity appear "staggered", i.e., odd-even os- 
cillations are visible, reflecting the finite excitation gap 
at unitarity. In the non- interacting limit, in contrast, the 
excitation gap is zero and no odd-even staggering of the 
hyperradial potential curves is visible. 

To extrapolate to the large A^ limit. Fig. [5] shows the 
normalized coefficients Cn, Eq. ([5^ with E^j replaced 
by Eni^etf, as a function of A^ (just as in our analysis 
of the energies i^oo we find it useful to smooth the en- 
ergies Enj). The coefficients Cn oscillate between two 
smooth curves, a curve for even N (circles) and a curve 



14 



0.55 

0.5 
0.45 

0.4 
0.35 

0.3 
0.25 

0.2 



xxxxxxxxx^x 



0.55 
0.5 
" 0.45 
0.4 
0.35 



><So<XxXX XX X 



0.05 0.1 

l/N 



0.15 







10 



15 20 25 30 

N 



FIG. 9: (Color Online) Normalized coefficients Cjv, Eq. (|32|) 
with Eni replaced by Eni,etf, as a function of A'^; values for 
even A'^ are shown by circles and values for odd by crosses. 
The dash-dotted line shows the value ^ — 0.42 obtained by 
FN-DMC calculations for the homogeneous system [stI . [gJ, 
while a dashed curve shows the value ^ — 0.508 obtained with 
a renormalization procedure [tJ- The inset shows the same 
quantities as a function of instead of A'^. 



for odd N (crosses). As A'' increases, the difference be- 
tween the two curves decreases. In the large A^-limit, 
the value of Cn for two-component Fermi gases at uni- 
tarity should approach the universal parameter ^ [20j . 
This can be shown by relating the ground state energy 
obtained within the hyperspherical framework, Eq. (P7)) . 
to the LDA prediction (see above), or by applying renor- 
malized zero-range interactions within the hyperspheri- 
cal framework [53 • The dash-dotted and dashed lines in 
Fig. [HI show the ^ value obtained by FN-DMC calcula- 
tions for the homogeneous system (^ = 0.42) [s^, H^l and 
the ^ value obtained with a renormalization procedure 
(^ = 0.508) [t^], respectively. It is generally believed 
that the FN-DMC calculations provide the most reliable 
estimate for ^ to date. For comparison, our energies for 
the trapped system predict £,tr = 0.467 (see above). The 
circles in Fig. [5] approach this value. We attribute the 
fact that £,tr is larger than the corresponding value of the 
bulk system, i.e., ^ = 0.42, to the comparatively small 
system sizes (A^ < 30) included in our analysis. If this 
was true, we would expect the circles in the main part of 
Fig. [5] to turn around at larger A^ values. We note that 
we cannot rule out that the nodal surface entering our 
FN-DMC calculations might not be optimal. 

In addition to equal-mass unitary systems, we study 
small systems with unequal masses at unitarity. Fig- 
ure [TO] shows the excitation gap A(A^) for A = 3 at uni- 
tarity as a function of the mass ratio k. A (3) decreases 
from about O.Shuj for k = 1 to about 0.3huj for k — 8. 
A decrease of the excitation gap as a function of k has 
recently also been reported for the homogeneous unequal- 
mass system at unitary [t^. To better understand the 
decrease of A(A^) with increasing k, triangles and squares 




FIG. 10: (Color Online) Circles show the excitation gap A(A'^) 
for A'^ = 3 as a function of the mass ratio k at unitarity. 
Triangles and squares show the chemical potentials ni (3) and 
/i2(3), respectively. 



in Fig. [TO] show the chemical potentials /ii(3) and /i2(3) 
for the two species. The decrease of fii is related to 
the fact that trimers with negative energy form for suf- 
ficiently large k. We additionally note that the densities 
of the light and heavy particles do not fully overlap. This 
effect is unique to the trapped system (the study of the 
homogeneous system with unequal masses [zll assumes 
equal densities of the two components and full pairing). 
Simple arguments lead one to conclude that a partial 
density overlap as opposed to a full density overlap leads 
to a decrease of the excitation gap. Thus, it is not clear 
if the decrease of A(3) visible in Fig. [TO] with k is due 
to the same mechanisms that lead to a decrease of A in 
the homogeneous system or due to the specifics of the 
trapping potentials, or possibly both. 



C. Excitation spectrum at unitarity 

Excitation spectra of two-component Fermi gases are 
rich. For four equal-mass fermions, e.g., Ref. [59j shows 
how the spectrum evolves from the non-interacting limit 
for small \as\, as < 0, to different families in the small 
as region, as > 0: One family consists of states that 
describe two bound dimers, another consists of states 
that describe a bound dimer plus two atoms, and yet 
another consists of states that describe a gas. Between 
these two limiting cases is the unitary region where the 
eigenspectrum is expected to be characterized by unique 
properties, similar to those of the non-interacting system 
(see Sec. Ill C"|) . In particular, in the unitary regime fam- 
ilies of eigenenergies separated by 2hu! are expected to 
exist [131. This prediction has recently been verified for 
up to six particles with equal masses to within numeri- 
cal accuracy, i.e., to within 2% [2^. Here we extend our 
analysis to unequal-mass systems with A^ = 4 and L = 0. 

Circles in Fig. [11] show the zero angular-momentum 
energy spectrum calculated by the CG approach for four 
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particles at unitarity as a function of k. The range of 
the Gaussian potential is Rg = O.Ola^Q. To analyze the 
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TABLE IV: CoefRcients of the hyperradial potential curves 
Vs^ (R), Eq. (1211, for the TV = 4 system with L = for various 
mass ratios k. 
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2.03 


4.46 


5.05 


8 


2.45 


3.81 


5.29 


2 


2.09 


4.41 


4.88 


9 


2.45 


3.74 


5.35 


3 


2.18 


4.27 


4.90 


10 


2.42 


3.68 


5.39 
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2.27 


4.15 


4.98 


11 


2.37 


3.62 


5.39 
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2.34 


4.04 


5.06 


12 


2.29 


3.57 


5.30 
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2.40 


3.95 


5.15 
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FIG. 11; (Color Online) Four-body energy spectrum for L — 
at unitarity as a function of k. Circles correspond to the 
numerical results obtained by the CG approach. Solid, dashed 
and dash-dotted lines show the energies E,jQ + 2nhuj for v = Q, 
1 and 2, respectively (n = 0, 1, • ■ ■ ). 

eigenenergies, we employ the hyperspherical framework. 
Assuming that the separation of the wave function (see 
Sec. Ill CI) holds for the short-range interactions consid- 
ered here, we expect that the energy spectrum consists 
of families of energy levels separated by 2fiw. Solid lines 
show the energies i?oo + 2ri?iaj [n non- negative integer), 
where i?oo denotes the lowest positive energy of the spec- 
trum (for sufficiently large k, negative energy states form; 
these are not shown in Fig. [TT|) . The agreement between 
the solid lines and the CG energies indicates that the 
2tiuj spacing, predicted for zero-range interactions, is ful- 
filled within our numerical accuracy. We repeat this ex- 
ercise for the next family of energy levels: Dashed lines 
show the energy Ei^ -\- 2nfujj, where Eiq corresponds to 
the lowest positive energy not yet assigned to a family. 
Similarly, dash-dotted lines connect states belonging to 
the third family. In addition to the just outlined assign- 
ment of quantum numbers, we checked in a few cases that 
the hyperradial wave functions {R) corresponding to 
the energies E^n possess n hyperradial nodes (see also 
Sec. lIVE| . The lines in Fig. [TT] show a crossing of energy 
levels belonging to different families at k w 4. In close 
vicinity to this crossing, the spacing may not be exactly 
2hjj. 

As already pointed out in the previous section, the en- 
ergies EijQ determine the coefficients Sy of the hyperradial 
potential curves Vs^{R). Table [TVl summarizes the three 
smallest coefficients for various n. To the best of our 
knowledge, these are the first calculations of the for 
four-particle systems with unequal masses. 

D. Structural properties along the BEC-BCS 
crossover 

In addition to the energetics, we analyze the one- 
body densities and pair distribution functions of two- 
component Fermi systems in the crossover regime for dif- 



ferent K. While the densities Pi{r) oi L — Q states are 
spherically symmetric, those of states with L > are not 
spherically symmetric. In the following, we determine 
the averaged radial densities Pi{r), normalized so that 
47r / pi{r)r'^dr = Ni] inr^ pi{r) /Ni tells one the proba- 
bility of finding a particle with mass at a distance r 
from the center of the trap. If A^i = A'2 and toi = m2, 
the radial one-body densities pi{r) and P2{r) coincide. If 
TOi and m2 or TVi and N2 differ, however, the radial one- 
body densities ^1(7") and P2i'f) are, in general, different. 
We also determine the averaged radial pair distribution 
functions Pij{r), normalized so that 47r J Pij{r)r'^dr = 1; 
4'Kr^Pij{r) tells one the probability to find a particle of 
mass jTii and a particle of mass mj at a distance r from 
each other. For notational simplicity, we refer to the 
radial one-body densities as one-body densities and to 
the radial pair distribution functions as pair distribution 
functions in the following. 

Figure [T^ shows the pair distribution function ^12(7') 
for iV = 3 (dash and dash-dotted lines correspond to 
L — and 1, respectively) and iV = 4 (solid lines) along 
the crossover for k = 1. Panel (a) shows results for 
as = — a/io, panel (b) for l/ag = and panel (c) for 
Us = O.la/io- Interestingly, the pair distribution func- 
tions for iV = 3 and 4 show a similar overall behavior. 
In the BCS regime [Fig. [T^ a)]. the quantity Pi2{r)r^ 
shows a minimum at small r (for very small r, Pi2{r)r^ 
goes smoothly but steeply to zero; this rapid change of 
Pi2(7')r^ is hardly visible on the scale shown in Fig. [T^ . 
At unitarity [Fig. [Hfb)], Pi2{r)r^ shows a maximum at 
small r and a second peak at larger r. In the BEC regime 
[Fig. [T2{c)], the two-peak structure is notably more pro- 
nounced. The peak at small r indicates the formation 
of tightly-bound dimers (one dimer for A'^ = 3 and two 
dimers for A = 4), while the peak between laho and 
2aho is related to the presence of larger atom-atom dis- 
tances set approximately by the atom-dimer distance for 
the three-body system and the dimer-dimer distance for 
the four-body system. This interpretation suggests that 
the three-particle system has one small and one large 
interspecies distance, and the four-particle system has 
two small and two large interspecies distances. Indeed, 
integrating Pi2{r) for A" = 3 and 4 from to the r 
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FIG. 12: (Color Online) Pair distribution functions Pi2{r), 
multiplied by r'^, for equal mass two-component Fermi sys- 
tems with TV = 3 and L = (dashed lines), A'^ = 3 and 
L = 1 (dash-dotted lines), and A'^ = 4 and L = (solid 
lines) obtained by the CG approach for three different scat- 
tering lengths ttsi (a) as = — flho (BCS regime), (b) l/us =0 
(unitarity), and (c) = O.laho (BEG regime). The pair dis- 
tribution function for A'' = 4 and as — O.laho [solid line in 
panel (c)] is shown in more detail in Fig. 1131 



value at which Pi2(r)r'^ exhibits the minimum, we find 
that the likelihood of being at small distances (forming 
a molecule) and being at large distances is the same. 

We now analyze the pair distribution function Pi2{t') 
for = 4 more quantitatively. Dash-dotted lines in 
Figs. fTSTa) and (b) show the pair distribution function 
Pi2(r), multiplied by r^, for two trapped atoms with — 
O.laho (normalized to 1/2). This dimer curve is essen- 
tially indistinguishable from the small r part of the four- 
particle pair distribution function (circles). To describe 
the large r part of the four-particle pair distribution func- 
tion, we consider two bosonic molecules of mass 2m, 
which interact through an effective repulsive potential 
with dimer-dimer scattering length add ~ 0.6as [2ll[23j. 
The dashed line in Fig. [THf a) shows the pair distribu- 
tion function for this system under external confinement. 
This dashed curve is essentially indistinguishable from 
the large r part of the pair distribution function for the 
four-particle system. For comparison, a dotted line shows 
the pair distribution function for two non-interacting 
trapped bosons of mass 2m. Figure [13] indicates that 
the effective repulsive interaction between the two dimers 
is crucial for reproducing the structural properties of the 
four-body system accurately. Our analysis shows that the 
entire pair distribution function P12 (r) of the four-body 
system in the weakly-interacting molecular BEG regime 
can be described quantitatively in terms of a "dimer pic- 
ture" . 

We now return to Fig. [12] and discuss how the 
symmetry-inversion of the A^ = 3 system along the 
crossover (see Sec. IIV A[) is reflected in Pi2{r). In the 
BCS regime and at unitarity [Figs.fTWa) and (b)], Pi2{r) 
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FIG. 13: (Color Online) (a) Circles show the pair distri- 
bution function Pi2(r), multiplied by r'^, for as — O.laho 
(BEG regime) calculated by the GG approach for A'^ = 4 and 
n = 1 [note, this quantity is also shown by a solid line in 
Fig. (T^Jc)]. For comparison, the dash-dotted line (blue on- 
line) shows Pi2(r)r'^ for two atoms of mass m with the same 
scattering length but normalized to 1/2, the dashed line (red 
online) shows Pi2(r)r^ for two trapped bosonic molecules of 
mass 2m interacting through a repulsive effective potential 
with add ~ 0.6as, and the dotted line (green online) shows 
Pi2{r)r^ for two trapped non- interacting bosonic molecules 
of mass 2m. Panel (b) shows a blow-up of the small r region. 



shows less structure for L = 1 than for L = 0. In the 
weakly-interacting molecular BEC regime [Fig. [T^ c)]. 
the pair distribution function for L = nearly coincides 
with that for L = 1 at small r but is more compact than 
that for L = 1 at large r. 

Next, we analyze how the behaviors of the pair distri- 
bution functions Pi2{r) ioi N = 3 and 4 change along the 
crossover if the mass ratio is changed from k = 1 to 4. 
Figure [H] shows the pair distribution functions for k = 4. 
For N = 3, we consider three-particle systems with either 
a spare light particle or with a spare heavy particle. The 
pair distributions for the three-particle system with two 
light particles and one heavy particle are notably broader 
than those for the three-particle system with one light 
particle and two heavy particles. This behavior can be 
attributed to the fact that aJ^^J > a^i^^. Besides this, a 
comparison of the pair distribution functions shown in 
Fig. [M] for K = A and those shown in Fig. [12] for k = 1 
reveals that the overall behavior of the Pi2{r) is similar. 

Figures [T5f a) and (b) show the one-body densities for 
K = 1 and 4, respectively. In the non-interacting limit 
[the solid lines show pi{r) and the circles show P2 (?")], 



the sizes of pi{r) and P2(f) are determined by a^/^^ and 
(2) 

aj-J , respectively. As is evident in Fig. [151 the density 
of the light particles extends to larger r than the den- 
sity of the heavy particles. The density mismatch for 
K = 4 between the two one-body densities decreases as 
as is tuned through the strongly-interacting regime to the 
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FIG. 14: (Color Online) Pair distribution function Pi2(r), 
multiplied by r^, for two-component Fermi gases with k — 
4 for different scattering lengths as: (a) as — —aho (BCS 
regime), (b) 1/as = (unitarity), and (c) as = O.laho (BEC 
regime). Dashed and dash-dotted lines show Pi2(r)r^ for TV = 
3 (two heavy particles) with L = and 1, respectively. Circles 
and squares show P\_2{i')r^ for N = 3 (two light particles) with 
L — and 1, respectively. Solid lines show Pi2{r)r^ for N = -i 
with L = 0. 

weakly-interacting molecular BEC side. In the weakly- 
interacting molecular BEC regime, two molecules con- 
sisting each of a heavy and a light particle form. In this 
regime, the size of the system is determined by the molec- 
ular trap length and the densities pi (r) and p2 (r) [trian- 
gles and dash-dotted line in Fig. [TSTb)] nearly coincide. 
Furthermore, the densities are to a very good approxima- 
tion described by the one-body density for two bosonic 
molecules of mass mi -I- m2 interacting through an ef- 
fective repulsive interaction characterized by the dimer- 
dimer scattering length {add ~ 0.77as for k = 4 pll[23|). 

We have also analyzed the pair distribution functions 
Pii{r) for K = 1 and 4 (not shown). The small r region 
of the Pii{r) is controlled by the Pauli exclusion princi- 
ple between identical fermions. In the weakly-interacting 
molecular BEC regime, the pair distribution functions 
Pii(r) and P22{f) nearly coincide even for k = 4. In this 
regime, the pair distribution functions Pair) are well ap- 
proximated by that for two particles of mass mi -I- m2 
interacting with a repulsive potential characterized by 

a-dd- 



E. Structural properties at unitarity 

This section discusses selected structural properties of 
two-component equal mass Fermi gases at unitarity with 
up to = 20 atoms. For small systems [N < 6), 
we present structural properties calculated using both 
the CG and the FN-DMC methods. For larger systems, 
however, our interpretation relies solely on the structural 
properties calculated by the FN-DMC method. 
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FIG. 15: (Color Online) One-body densities pi(r) and P2{r) 
for A'^ = 4 and (a) k = 1 and (b) k = 4 for different scattering 
lengths as [for k, = 1, pi(r) and p2(r) coincide and only p2(r-) 
is shown]: Circles and solid lines show pi(r) and p2(r) for 
as — 0, squares and dashed lines show pi(r) and P2(r) for 
1/as =0, and triangles and dash-dotted lines show pi(r) and 
P2(r) for as = Cla^o. Note, p2{r) for k = 4 and as = [solid 
line in panel (b)] is multiplied by a factor of three to enhance 
the visibility. 

To assess the accuracy of the nodal surfaces employed 
in our FN-DMC calculations as well as of the accu- 
racy of the mixed estimator [see Eq. (|4T|) in Sec. IIII Bj . 
Figs. llGf a) and (b) compare the pair distribution func- 
tions Pi2{r) for the three-particle system with L — 1 
and the four-particle system with L = 0, respectively, 
calculated by the CG and the FN-DMC methods. The 



(a) 




15 3 4 5 



FIG. 16: (Color Online) Pair distribution functions Pi2(r), 
multiplied by r^, at unitarity for equal mass Fermi systems 
with (a) iV = 3 (L = 1) and (b) = 4 (L = 0) atoms 
calculated by the CG method (solid lines) and by the FN- 
DMC method (circles). The agreement is excellent. 

agreement between the pair distribution functions calcu- 
lated by the CG method (solid lines) and the FN-DMC 
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method (circles) is very good, validating the construction 
of the nodal surface of tpx- Furthermore, the good agree- 
ment suggests that the mixed estimator results, for the 
guiding functions employed, in structural properties very 
close to those one would obtain by an exact estimator. 

Figure [17] shows the pair distribution functions Pi2{r) 
calculated by the FN-DMC method for equal mass Fermi 
systems with TV = 3 — 20 at unitarity. To simplify the 




FIG. 17: (Color Online) Dashed and solid lines show the pair 
distribution functions Pwir), multiplied by r^Npair {Npair 
denotes the number of interspecies distances) , for equal mass 
Fermi systems at unitarity with even N (TV = 4, 6, ■ ■ ■ , 20) 
and odd TV (TV = 3, 5, ■ • • ,19), respectively, calculated by 
the FN-DMC. Beyond r ^ aho, Pi2{r)r^ Npair is smallest for 
TV = 3 and largest for TV = 20. 

comparison. Fig. [T71 shows the even TV results as a dashed 
line and the odd TV results as a solid line. Furthermore, 
Pi2 {r)r'^ is multiplied by the number Npair of interspecies 
distances so that the TV = 3 distribution function has the 
smallest and the TV = 20 distribution function the largest 
amplitude for r > a^o- The pair distribution functions 
for even TV show a similar behavior for all TV considered; 
both the small r and the large r peaks grow monoton- 
ically and smoothly with increasing TV. For odd TV, in 
contrast, the small r peak changes somewhat discontin- 
uously at TV « 11. This behavior can be attributed to 
the guiding functions employed. For even TV, the guiding 
function ipTi, whose nodal surface is constructed from 
the two-body solution, gives the lowest energy for all TV 
(except for TV = 4). For odd TV, however, ipT2 results 
in a lower energy for TV < 9, ?/'T3 for TV = 11, and tpxi 
for TV > 13. Thus, the pair distribution functions clearly 
reveal how the structural properties depend on the nodal 
surface employed in the FN-DMC calculations and pro- 
vide much deeper insights into the different ipx employed 
than a mere comparison of the energies. 

For TV > 13, Fig. [T7| indicates that the amplitudes 
of the scaled interspecies pair distribution functions are 
nearly the same for neighboring systems. For example, 
the quantities Pi2[r)r'^ Npair for TV = 18 and 19 agree 
to a good approximation, suggesting that one can think 



of the TV = 18 system as consisting of nine pairs, and 
of the TV = 19 system as consisting of nine pairs plus 
a spare atom. Note that this interpretation hinges crit- 
ically on the nodal surface employed in our FN-DMC 
calculations; a small change in the nodal structure of the 
guiding function may change the small r behavior of the 
pair distribution functions non-negligibly. 

We next investigate in Fig. [TSl where the spare particle 
is located in the odd- TV systems at unitarity. This fig- 
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FIG. 18: (Color Online) The one-body density pi(r) (solid 
lines) is shown for TV = 3, 9 and 15 (for r > 0.5aho, from bot- 
tom to top), together with the one-body density p2{r) (dashed 
line) for TV = 3, 9 and 15 (for r > O.^aho, from bottom to top) 
for equal mass two-component Fermi gases at unitarity calcu- 
lated by the FN-DMC method. 

ure shows the one-body densities pi[r) (solid lines) and 
P2{f') (dashed lines) for TV = 3, 9 and 15. For TV = 3, the 
difference between pi{r) and P2{f) is roughly constant 
across the trap. The behavior is similar for TV = 9. In- 
terestingly, the densities for TV = 9 show a minimum at 
r = 0, reflecting the fact that the FN-DMC calculations 
employ the nodal surface of the ideal Fermi gas, i.e., use 
'0T2 [Eq. (|40)) ]. For TV = 15, the nodal surface employed 
is constructed from the two-body solution [see Eq. ([55]) ]. 
and consequently, the behavior of the density profiles dif- 
fers from that for the smaller TV. The densities pi{r) and 
P2{f') nearly coincide at small r. At large r, however, the 
density pi{r) has a larger amplitude than p2{i') (recall 
TVi = TV2 -I- 1). Our data for TV = 15 indicate that the 
spare particle is not distributed uniformly throughout the 
trap but has an increased probability to be found near 
the edge of the cloud. Possible consequences of this find- 
ing for the excitation gap have already been discussed in 
Refs. HO, [73. 

To quantify the analysis of the one-body densities, we 
integrate piijr) over r, 



N,{r) = An / pi{r'y^dr' 



(44) 



For a finite upper integration limit, TVi(r) monitors how 
many of the TV^ particles are located between zero and 
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r. Figure fT9l shows Ni{r) for iV = 3, 9 and 15. As in 
Fig. [THl the results for component one are shown by sohd 
lines and those for component two by dashed lines; in the 
large r limit, the Ni{r) equal Ni^ as expected. One can 
now read off nicely, in which r-regions the densities of 
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FIG. 19: (Color Online) Solid and dashed lines show the in- 
tegrated quantities N\{r) and N2(r) [Eq. (|44p ]. respectively, 
as a function of r for a two-component Fermi gas at unitar- 
ity. At large r, the curves correspond from bottom to top to 
iV = 3, 9 and 15. 

the two components agree and where they disagree. For 
iV = 3, e.g., the two atoms of component one and the one 
atom of component two are added over approximately 
the same r-region. For N = 15, in contrast, the first five 
atoms of the two components are located in the region 
with r < 1.5aho', this core region can be considered "fully 
paired" . The last three atoms of component one and the 
last two atoms of component two form, loosely speaking, 
a "partially paired or unpaired outer shell" . We find 
similar behaviors for the odd-A^ systems with A^ = 13, 17 
and 19. It will be interesting to see if this interpretation 
holds for larger A^, and if this information can be used 
to shed light on the phase diagram of asymmetric Fermi 
gases HMI- 

To further verify the validity of the special proper- 
ties of two-component Fermi gases at unitarity as well 
as to further assess the accuracy of our guiding func- 
tions employed in the FN-DMC calculations, we analyze 
the hyperradial densities ^00(2;) for various A^. Symbols 
in Fig. [20] show the dimensionless hyperradial density 
Fqq{x) calculated using the mixed Monte Carlo estima- 
tor, Eq. (I4ip . for A^ = 3 to 10. Here, x is the dimension- 
less hyperradius defined just above Eq. l|30|) and the nor- 
malization of the ^00 is chosen so that FQQ{x)dx — 
1. The dimensionless hyperradius x is scaled by R'pfj, 
which we evaluate by approximating Epfj in Eq. ([?T|) by 
En I, EFT- This is similar to the "smoothing procedure" 
discussed in the context of Figs. [7] and O The hyperradial 
densities become more compact as A^ increases, owing to 
the increase of the effective mass fief / entering into the 
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FIG. 20: (Color Online) The hyperradial density F^oi^) is 
shown as a function of the dimensionless hyperradius x for 
A = 3 — 10, calculated using the mixed Monte Carlo esti- 
mator (symbols) and the analytical expression with the FN- 
DMC energies (lines), respectively. The maximum of Fqq{x) 
is smallest for A = 3 and largest for N = 10; the MC results 
are shown by filled circles for A = 3, open circles for A = 4, 
filled squares for N — 5, open squares for A = 6, filled trian- 
gles for A = 7, open triangles for N — 8, filled diamonds for 
A = 9 and open diamonds for A = 10. 



effective hyperradial Schrodinger equation [Eq. ([50)1 ] with 
increasing A^. Furthermore, the maximum of the hyper- 
radial densities occurs at slightly larger x values for odd 
A^ systems than for even A^ systems, in agreement with 
the odd-even staggering discussed in Sec. lIVBl in the con- 
text of the hyperradial potential curves V{R). 

In the limit of zero-range interactions, the adiabatic 
approximation is expected to be exact (see Sec. Ill C|) . In 
this case, the functional form of the hyperradial wave 
functions is known analytically [see Eq. (^5]) ]. and can 
be compared with the Monte Carlo results obtained for 
short-range potentials by sampling the total wave func- 
tion and integrating over all coordinates but the hyperra- 
dius. Solid lines in Fig. [^D] show the hyperradial densities 
Fqq{x) for A^ = 3 to 10 predicted analytically for zero- 
range interactions, using the FN-DMC energies Em listed 
in Table [m] of this paper and Table II of Ref. The 
agreement between the analytical results and the Monte 
Carlo results obtained for finite range potentials is quite 
good. On the one hand, this agreement lends numerical 
support for the separability or near separability of the to- 
tal wave function into a hyperradial and a hyperangular 
part. On the other hand, the good agreement suggests 
that the nodal surface employed in our MC calculations 
is appropriate. 

Finally, we analyze the hyperradial densities calculated 
by the CG approach for the N = 6 system. Since the CG 
approach allows for the determination of excited states, 
this analysis allows us to verify that the 2huj spacing re- 
ported in Ref. [23| and discussed in Sec. Ill CI corresponds 
indeed to breathing-mode excitations, i.e., to excitations 
along the hyperradial coordinate. To extract the hyper- 
radial densities, we integrate the square of the wavefunc- 
tion 5''''^' over all the coordinates but the hyperradius. If 
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FIG. 21: (Color Online) Hyperradial density F§„{R) for n = 
1,2 and 3. Here, we choose ^jv = m {£, = auo)- The solid 
lines show the analytical solutions while the circles show the 
numerical results obtained by integrating {^'^"^Y' calculated 
by the CG method over all coordinates but the hyperradius 
R. 

the universal behavior is fulfilled, then the hyperradial 
densities should coincide with the square of the analyt- 
ically determined F^n{R), Eq. (|28|) . which are shown in 
Fig- EH by solid lines. The integration over the hyperan- 
gular coordinates is carried out using Monte Carlo inte- 
gration techniques. Symbols in Fig. 1211 show the result- 
ing hyperradial densities for ^ = and n = 0, 1 
and 2. The numerically determined hyperradial densi- 
ties indicate that the excitations are to a good approx- 
imation located along the R coordinate, supporting the 
interpretation of the 2hu} spacing within the hyperspher- 
ical framework. The agreement between the numerical 
and analytical results is excellent for the ground state 
[see Fig. [STTa)]. For the excited states with n = 1 and 2, 
the small deviations between the numerical and analyti- 
cal results may be due to finite range effects or not fully 
converged numerical results. 

V. CONCLUSIONS 

This paper presents a microscopic picture of the prop- 
erties of ultracold two-component fermionic systems in a 
trap. Complementing previous studies [13, |2l[, we focus 
on the energetics of odd N systems, and the structural 
properties of both odd and even N systems. 

For sufficiently few particles, we solve the Schrodinger 
equation for equal and unequal mass systems, starting 
from a model Hamiltonian with short-range interspecies 
s-wave interactions using the CG approach. This ba- 
sis set expansion technique allows for the determination 
of the eigenspectrum and eigenstates with controlled ac- 
curacy throughout the BEC-BCS crossover. We find 
that the spectrum and the structural properties of small 
trapped two-component Fermi systems change qualita- 



tively throughout the crossover regime. 

An analysis of the energies of the = 3 systems in 
the weakly-interacting BEC and BCS regimes allows us 
to determine the validity regime of the analytically de- 
termined perturbative expressions for small |as|. Fur- 
thermore, we find that the angular momentum of the 
iV = 3 ground state changes from L — 1 in the weakly- 
attractive BCS regime to L = in the weakly-repulsive 
BEC regime for all mass ratios considered. By addition- 
ally treating the N — 2 and 4 systems, we determine 
the excitation gap A (3) throughout the crossover region: 
For equal frequencies, the excitation gap decreases for all 
scattering lengths with increasing mass ratio. For TV = 4 
systems with k < 13, we determine the L — excitation 
spectrum at unitarity. The spectrum determines the s^, 
coefficients of the hyperradial potential curves and also 
verifies within our numerical accuracy the 2^01 spacing 
prediction, which was derived analytically assuming uni- 
versality [5J|. We verified in a number of cases that the 
2hLL! spacing corresponds indeed to breathing-mode exci- 
tations. 

Our analysis of the energetics is complemented by 
studies of the structural properties. For the four-particle 
system with equal and unequal masses, e.g., we show how 
the pair distribution functions in the Os > region (small 
tts) can be described by a system of two molecules inter- 
acting through an effective dimer-dimer potential with 
positive dimer-dimer scattering length. A similar analy- 
sis was carried out for the A^ = 3 system and we verified 
that this system behaves as an interacting system of an 
atom and a dimer. 

Our small A^ studies have implications for optical lat- 
tice experiments. Our results can be applied directly 
if each optical lattice site is approximately harmonic in 
the non-tunneling regime. In this context, Ref. [sij al- 
ready pointed out, including the energies of the two- and 
three-particle system, that the occupation of optical lat- 
tice sites with three equal-mass atoms is unfavorable. To 
start with let us consider a system with two lattice sites 
and six atoms (three of each species). We imagine that 
the lattice sites are loaded by adiabatically turning up 
the barrier height between the two sites. It follows from 
our energies calculated for A^ = 2 through 4, that the 
system ends up with unequally occupied lattice sites at 
the end of the ramp: For both equal and unequal masses, 
and for all scattering lengths, the energy of the two-site 
lattice is minimal if two atoms occupy one site and the 
other four atoms occupy the other site. This is simply a 
consequence of the fact that A (3) is positive throughout 
the crossover (note, the energy of the unequal mass sys- 
tem might be further lowered if we consider the formation 
of pentamers and sextamers with negative energies; these 
states are not included in our analysis). The arguments 
presented here for just two lattice sites generalize readily 
to lattices with multiple sites. 

Instead of ramping up the barrier height adiabatically, 
we now imagine a fast non-adiabatic ramp. In this case, 
the likelihood of finding three atoms per lattice site at 
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the end of the ramp is finite. Since the excitation fre- 
quencies for two-, three- and four-particle systems are 
different, a "purification sweep" [lB| can be used to then 
prepare a system with either three or no particles per 
site. These three-particle systems could be investigated 
spectroscopically (see, e.g.. Sec. IIV Cl for a discussion of 
the excitation spectrum of the four-particle system) . Al- 
ternatively, one might ask whether it would be possible 
to measure the odd-even physics by adiabatically lower- 
ing the lattice barrier and monitoring the point at which 
tunneling sets in. 

In addition to systems with equal numbers of atoms in 
the two species, we consider an optical lattice with twice 
as many heavy as light atoms. If the mass ratio is suffi- 
ciently large, trimers consisting of two heavy atoms and 
one light atom with negative energy can form at each 
lattice site, paving the way for spectroscopic studies of 
these delicate systems. Furthermore, by starting with a 
bound trimer in a deep lattice and then lowering the lat- 
tice height, a gas consisting of bound trimers can possibly 
be prepared. 

To extend the studies of the energetics and structural 
properties to larger systems, we employed the FN-DMC 
technique. This approach determines the lowest energy of 
a state that has the same symmetry as a so-called guiding 
function and thus an upper bound to the exact eigenen- 
ergy. Detailed comparisons of the energies and the struc- 
tural properties calculated by the CG and FN-DMC ap- 
proaches benchmark the nodal surfaces employed for sys- 
tems with N < 6. In the strongly-correlated unitary 
regime, e.g., the FN-DMC energies for equal- mass two- 
component Fermi systems agree with the CG energies to 
within 2%. 

Our even N energies {N < 30) for equal-mass sys- 
tems at unitarity show vanishingly small shell structure. 
Applying the local density approximation and approxi- 
mating the non-interacting energies by the corresponding 
extended Thomas-Fermi expression, we find ^tr — 0.467, 
which is somewhat larger than the value of the homo- 



geneous system, ^^om = 0.42. We note that the expres- 
sion ^/^Eni^eft describes the equal-mass energies for 
iV < 30 at unitary very well; the small disagreemet be- 
tween ^tr and S,hom is most likely due to the small number 
of particles considered in the present work. Combining 
the even and odd N energies, we find that the excita- 
tion gap A(Ar) at unitarity increases with N. Also, the 
one-body densities and pair distribution functions at uni- 
tarity are studied for up to iV — 20. For odd N with 
N > 11, we observe that the extra "unpaired" particle 
is located predominantly near the edge of the cloud, in 
agreement with previous predictions pol [TOj . This sug- 
gests that the LDA cannot be applied to determine the 
excitation gap. Furthermore, we find that the hyperra- 
dial densities of the lowest gas-like state with iV < 10 
calculated for short-range interactions by the FN-DMC 
method agree with the analytically predicted ones, indi- 
cating that the lowest gas-like state does indeed behave 
universally. Selected hyperradial densities for larger N 
were already presented in Ref. [20j . 

The energies and structural properties for the equal- 
mass two-component Fermi systems at unitarity pre- 
sented in this paper may serve as a benchmark for other 
calculations. Recently, e.g., a DFT treatment determined 
the energies for systems with up to iV = 20 particles [ll] . 
The good agreement between the FN-DMC energies and 
the DFT energies suggests that the functional employed 
in the DFT calculations captures the key physics. How- 
ever, close inspection of the FN-DMC and DFT energies 
indicates that the agreement between the even N and 
odd N energies is not equally good. While this could be 
a consequence of the nodal surfaces employed in our FN- 
DMC calculations, it could alternatively indicate that the 
DFT treatment employed in Ref. [zl] for odd N is not 
optimal. Thus, it is hoped that our results will proof 
helpful in assessing the accuracy of the DFT approach 
and other approaches. 

We acknowledge support by the NSF, and fruitful dis- 
cussions with J. D'Incao and S. Giorgini. 
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